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Abstract Simulations of motion by mean curvature in bounded domains, with applications to bubble motion 
and grain growth, rely upon boundary conditions that are only approximately compatible with the equation of 
motion. Three closed form solutions for the problem exist, governing translation, rotation and expansion of a 
single interface [1], providing the only benchmarks for algorithm verification. We derive new identities for the 
translation solution. Then we estimate the accuracy of a straightforward algorithm to recover the analytical 
solution for different values of the velocity V given along the boundary. As expected, for large V the error can 
reach unacceptable levels especially near the boundary. We discuss factors influencing the accuracy and propose 
a simple modification of the algorithm which improves the computational accuracy. 

1 Introduction 

Foams are ubiquitous in everyday life [2]. They are used daily in the home, in both food and cleaning products. 
Moreover, their industrial uses are varied and contribute significantly to modern industrial processes. For 
example, foams are used to separate ores (e.g. zinc, lead) from the rock in which they are found, and to push 
oil out of porous rock. They are also used in the decontamination cleaning of vessels, and in firefighting. In 
all these applications, it is the flow of foam that is driving the process. Therefore, if it were possible to predict 
how a foam would behave in this range of scenarios, when subjected to a complex collection of deformations, 
each process could be made more efficient, either in terms of yield/output or cost-effectiveness. 

Aqueous foams and concentrated emulsions, which are very similar to foams, have peculiar and remarkable 
properties: they are complex fluids whose properties lie between the familiar extremes of liquid and solid. For 
small strains, a foam behaves as an elastic solid while at large strains (or strain-rate) a foam moves like a liquid. 
They therefore generate a rich range of behaviours, but with a well-defined local structure that obeys Plateau's 
geometric laws and the Laplace- Young law that balances interface curvatures with bubble pressure differences. 
Our idea is to use the precise, known, structure of a foam to predict its rheological response. 

One of the leading tools for the analysis of foam structure is the Surface Evolver developed by Ken Brakke 
[3]. This consists of software expressly designed for the modelling of soap bubbles, foams, and other liquid 
surfaces shaped by minimizing energy (such as surface tension and gravity), and subject to various constraints 
(such as bubble volumes). Originally designed to model grain growth in metals, it is easily extended to the case 
of foams, and has been used in a number of engineering disciplines to look at, for example, capillary surfaces, 
solder joints and fluid behaviour in microgravity. A surface is represented as a collection of triangles, so that 
the complicated topologies found in foams are routinely handled. In particular, the Evolver can deal with the 
topological changes encountered during foam flow. 

Here, we concentrate on a two-dimensional model of a foam, for case of analysis. We use the following 
"viscous froth model" (VFM) [4] to examine the evolution of foam films under shear: 

Pi -Pj =1 Kij + XVij (1) 

where 7 is the surface tension in the films (assumed constant), is the curvature of the film separating 
bubble i from bubble j and Vij its normal velocity. Each bubble has a well-defined pressure pi . In essence, the 
model augments the (equilibrium) Laplace- Young law with a term proportional to the velocity; the constant 
of proportionality A is a drag coefficient, representing the external dissipation due to friction with the walls of 



the container, allowing the investigation of strain-rate effects. If the external dissipation is negligible, then we 
recover the Laplace- Young law for a foam in equilibrium and a model for quasi-static evolution. 

A Surface Evolver simulation with the VFM proceeds in the following way. Each interface is discretised 
into short straight segments and curvatures K^j are calculated pointwise at the intersections of these segments. 
Bubble pressures are obtained by deriving a matrix equation based upon an integral of ([1]) around each cell. 
Then the endpoints of the segments are moved according to the motion equation ((T|) using a small time step 
At. Boundary conditions are applied according to the problem under consideration. 

In the ideal model of the evolution of crystalline grains in a polycrystalline metal, known as normal grain 
growth, the size of each grain evolves due to the normal motion of each of its boundaries [5] . Each boundary 
has a certain "mobility" A, and moves in such a way as to reduce the total perimeter of the pattern, as in foams. 
However, without area (volume in 3D) constraints, this is motion by mean curvature. As noted above, this is a 
well-posed limit of the VFM (see equation ([1])) when bubble pressure differences are negligible, such as in freely 
translating films (grain boundaries) and ordered (hexagonal) foams. The grain growth model requires that we 
have 120° at vertices, justifying a posteriori that assumption in the VFM. 

Here, we ask how such a solution can be commensurate with the boundary of the domain and how well this 
solution can be calculated numerically. To the best of our knowledge, answers to this question are not available 
in the literature. 

2 Curvature-driven motion of a bounded interface 
2.1 Problem formulation 

In vector form, the motion of an interface in the model of ideal grain growth ([T]) can be described by 



V 



KTlI, 



(2) 



where n and s are the normal and tangential unit vectors to the interface (see Fig. 1): 



n = [711,712], s [n2, -7ii]. 



(3) 



If the representation of the interface is taken in the form: 



x^x(y,t), y e[y{t),y{t)] 



(4) 



then the vector components tii, 712 are calculated as follows 




(5) 



where Xy = dx/dy, ds = y/ (dx)"^ + {dy^ and is the tangential angle to the interface (see Fig. 1). 
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Figure 1: The bounded interface considered here. 



Finally, the vector v = [wi,W2] is the instantaneous velocity of the point (x,y) lying on the interface at time t 
and K is the curvature of the interface at that point: 



Equation ^ can be also written in component form: 



w„ = V • n = viui + V2n2 = k, (7) 

Ws = V • s = — wirt2 + i'2'^i = 0. (8) 

In this paper we will consider only MuUins' solution [l], also known as the grim reaper, (see below) which 
is symmetrical with respect to the x-axis. Taking into account the direction of the interface motion, we can 
assume in what follows that: 

ni < 0, 712 > 0, Xy > 0, < 6 < n/2, Xyy > 0, k < 0, (9) 

Vn < 0, {V2 < 0, VI > 0). (10) 

Equation ([7]) is widely discussed in the literature [U |6], while the second equation ^ is somehow usually 
forgotten in this context. If one is only interested in reconstruction of the interface position at any time step 
an approach based only on equation ([7]) is sufficient. However, if it is required to control the position of each 
tessellation point along the interface, as in the case of numerical computation, then both equations are equally 
important. Note that there has previously been an attempt to control both the velocity components in a specific 
way in [5] . Equation ([H]) allows us to find relation between the two unknown components of the velocity vector 
V and the normal vector n in the form: 

(11) 



V2 

712 — — ni. 



This allows us to eliminate components of the normal vector n from equation ([7]) to give: 
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(12) 



2.2 Mullins' solution for translation revisited 

Let us assume that the interface conserves its shape but moves in the x-dircction with a constant speed V. We 
consider two points A and C having the same y-coordinate y = yo sX two consecutive time steps and to + dt 
(see Fig. 2). It is clear that these two points correspond to two different material points. Namely, there exists 
a point B on the interface at time to which moves according to the curvature law ^ to the point C for an 
infinitesimally small time step dt. If the coordinates of the point A are (a;o,i/o) then the coordinates of B and 
C can be written (xo + dx, yo + dy) and (xq + Vdt, yo). 



y, + dy 




+ dx 

Figure 2: Interface under translational motion at two consecutive instants in time to and to + dt. 
Note that 

BC = 'v{xo + dx, yo + dy)dt = \r{xo,yo)dt + 0{dtds). 
On the other hand, \BC\ = V sinOdt, tanO = vi/\v2\. As a result one can conclude: 
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{viiy) + vl{y)), 



(13) 



in some interval y € [0,h]. This relation immediately follows in the case of the translation movement of the 
interface in cc-direction from (jl2|) . 

Now, to reconstruct the solution obtained by MuUins [1] it is sufficient to substitute ^ into (fT^ to have: 

TT/2~Vy = e. (14) 

Here we have taken into account the second of the two symmetry conditions at the point y ~ 0: 

V2{0) = 0, 9{0)=n/2. (15) 

Equation fH]) can be written in the form 

cotVy = y^, (16) 
which after direct integration leads to MuUins' solution [1]: 

x(y) =x(0) + Ft- -^logcos(Fy), ye[0,h]. (17) 

where x{0) is the arbitrary initial position of the centre of the interface. This solution exists only under the 
condition h < hmax, where 

TT 

h,nax = (18) 

Note also that the angle 9 defined by such a solution monotonically decreases in the interval y €E (0, h) {h < 
hmax), taking values 

{h)^TT/2~Vh. (19) 
It is now possible to write analytical representations of all problem variables in the interval y € (0, h): 

Vi=Vcos^Vy, W2 = ~^l^sin2Fy, ni^smVy, n2—cosVy, K = VcosVy. (20) 

Note that the first symmetry condition (fT5|l has not been used but the reconstructed MuUins solution (fTT]) 
satisfies it automatically. It exhibits the following asymptotics near the symmetry axis: 

xiy)^x{0) + Vt+^Vy^ + Oiy''), y ^ 0; (21) 

thus near y = the interface is close to a parabola. Near the other end of the reaper (in the case of the maximal 
thickness h = hmax), the following asymptotic estimate can be obtained: 

x{y) ^ -^log{hmax " y) + 0(1), y ^ hmax. Or y - hmax ^ -de'^"" , a: ^ +00, (22) 
where d > is a constant. 

Remark: Note that relation (|13p also represents the boundary condition for any moving interface whose upper 
point lies on the line y = y = h and which moves in the .T-direction with velocity V . Moreover, the velocity can 
be in that case a function of time V = V{t): 

V{t)vi{h,t) = vl{h,t) + vl{h,t). (23) 

Substituting ([2]) into (|23p such a boundary condition can be equivalently rewritten in other forms: 

K{h,t) = V{t)ni{h,t), or Xyy{h,t)^V{t){l + {xy{h,t))^). (24) 

Note that we have not used pT]) to define ([24|. 



2.3 Arbitrary instantaneous solution of equations ([7j) - ([8]) in bounded domain. 

Let us consider any instantaneous solution of the equations ([Zl) - (|8]) with the prescribed boundary condition 
(|23p . Effectively this means that, for a particular time t, the end points of the interface y = y{t) and y = h 
are defined and the velocity components vi{y,t) and V2{y,t) are known functions of the variable y, while the 
problem is now to determine, using this information, the position of the interface in space variables {y,x). 
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Let us introduce a function which in what follows is considered known: 



Fiy)=vi{y) + 



>0, yeiy,h). 



(25) 



Note that the condition p3p may be not valid at all inside the interval y £ {y,h) as it was for MuUins' solution; 
as a result, F{y) is not a constant, in general. Equation (fT2l) can be integrated to give: 



(26) 



Here, recall that F depends upon time t, so that the constant of integration 9_ = 6_{t) and y = y{t). The equation 
(|26|) should be considered together with equation pT|) which, in this case, takes the form: 



dx 
dy 



V2 
Vl 



cot 0, 



or 



(27) 



(28) 



Equations (|26p and together indicate that the functions vi and V2 cannot be chosen arbitrarily to satisfy 
the vectorial equation ([2]) as one might expect. Instead, the following identity has to be satisfied: 



arctan— = / F{^)dC + e, 

V2 J 

y 

or writing w = Vi/v2, equation ([29]) becomes: 



(29) 



(30) 



which leads to the following identity valid within the entire interval y e (y, h): 



I 



vliy) = -w|(y) 



J 



(31) 



where the constant of integration clearly depends on time too, c = c(t). Note that any solution of equations ^ 
- ([8]) satisfies this additional relation, which makes sense only under the additional constraint: 



y 

1 < V2{Od^ + c<0. 



(32) 



As < 9 < 9, one can also deduce another condition which has to be true for any admissible velocities: 

y 



< J Fi^d^ < e. 



(33) 



Note that the constant x in (j28p is arbitrary (it changes only the position of the interface in the cc-direction 
and does not influence any other variables). To determine the other constants^(t) and c(t), we need to use the 
boundary conditions at the ends of the interface. Thus, condition (|23p together with (|3ip leads to 

(34) 



Vi{h) = V{l + C + 2l2), V2{h) = -VyJ- (1 + c + 2/2) (c + 2/2). 
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where we denote 

h 



h = h{t)= J V2iOdC. 
If the boundary condition on the other end is given in the form 

e{y{t))^m, (35) 

then all the constants have been defined. Such an instantaneous solution, assuming that the functions and 
V2{y) satisfy equation ([29|l . conditions ((34|) and restrictions p2|) and ([33]) . can be realised during the interface 
evolution at some step. 

In the case of the symmetrical solution, where both symmetry conditions psp and the additional condition 
vi{Q) — W > have to be satisfied, one can show that 

c(t) = 0, v2{y)^~W^y, y 0. (36) 

Note here that the value W = W{V, h) should be found from the constructed solution and is not an additional 
(arbitrary or given) constant. 

Finally, both restrictions (|32p and ([33|l should be valid for the symmetrical interface: 



n n 

J V2{m > -1/2, J F{Od^ < 7T/2. (37) 



Note that the tangential angle 9 for this solution is a monotonically decreasing function in the interval y € (0, h) 
so that 

h 

e{y) € {9,mn,TT/2), 0„„ = '"/^ " / ^(^^'^^^ ^^^^ 



It is straightforward to see that MuUins' solution ((20|) satisfies all these relationships with c{t) = 0, 9_{t) — 7r/2 
and W — V , as expected. In the Appendix, we construct analytical examples of symmetrical instantaneous 
solutions which are different from MuUins'. Those solutions are not, generally speaking, steady-state ones. 
This means that they can be reached at some time step i, given the interface boundary velocity V{t) and the 
position of the ends /i(t)), but all these parameters may later change with time. What is extremely interesting 
about these solutions is that some of them are well-defined for any velocity V > Q and an arbitrary position of 
the boundary y = h. This shows a "rich behaviour" of possible instantaneous solutions. It is also clear that 
there is an infinite number of admissible instantaneous solutions. Some of them can be realised during some 
specific non steady-state interface motion. For example, any instantaneous solution obtained during a numerical 
computation, for any particular time step, boundary velocity and topology, has to satisfy all the relations (|25p 
- ([35]). This will allow us to use the relations as indicators of the accuracy of computations. Moreover, they 
could provide a means to improve the accuracy of the algorithms. 

Note also that a family of arbitrary (asymmetric) instantaneous solutions constitutes an even larger set 
in comparison to the symmetric case. In fact, the family of symmetrical solutions has one degree of freedom 
(since the constant c{t) is equal to zero in this case) and correspondingly one less boundary condition (compare 
([35]) and (|15p ). Moreover, in the context of further applications of this result to a given algorithm, where the 
angle- type boundary condition has to be preserved at the interface intersection point, it is worth mentioning 
that the boundary condition ([35]) is therefore more important for application than the symmetry condition. On 
the other hand, symmetrical instantaneous solutions can also be considered a subset of the asymmetric solutions 
if one considers the interval {ho,h) instead of {0,h); {0 < ho < h). This idea has been exploited previously in 

3 Numerical Simulations 

To indicate the computational inaccuracy, we discuss MuUins' solution for a symmetrical reaper, for which all 
quantities arc known in closed form (see (jl4p and remarks thereafter), and compare it with the result of a 
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Figure 3: (a) The test problem considered liere consists of a single interface that is sheared symmetrically by translating 
the boundaries. The shape should correspond to Mullins' solution, (b) Standard algorithm for implementation of the 
boundary condition at each time step At. Aj and Bj are the respective tessellation points on the interface at times to 
and to + At. (c) Example of a multi-bubble foam simulation, for which the numerical procedure developed here will be 
of use. 

numerical computation using a simple algorithm (implemented in the Surface Evolver). This takes the form of 
a single interface separating two bubbles of equal pressure being pulled at a velocity V at each boundary (fig. 

131). 

The numerical procedure can be briefly described as follows. We start from a straight (vertical) line joining 
the two walls a distance 2h = 2 apart. This is subdivided into 2^ short elements ("edges" which meet at points). 
A time step At = 1 x 10~^ is chosen for the computations and bounds on the possible length L of each edge 
(0.01 < L < 0.05). The algorithm proceeds as follows: (i) each boundary point is moved in the x direction 
a distance VAt; (ii) the curvature at each point that is not on the boundary: k = F ■ n/L where L is the 
average length of the neighbouring edges and F is the negative energy (perimeter) gradient [3]; (iii) each point 
that is not on the boundary is then moved according to Ar = At k n (fig. [3]b)). This procedure (one "step") 
is repeated until the centre-point of the interface moves less than a critical value (1 x 10^^). Every 20 steps 
we check the edge length bounds and add or remove edges as necessary. Note that this standard algorithm 
preserves a reasonable restriction on the length of the edges; however, it works in a way that does not guarantee 
equal length edges. 

Note that this choice of the parameters for numerical simulation is standard and allows us to obtain accept- 
able accuracy in reasonable computational time [7j. On the other hand, when one computes the dynamics of 
foams with many bubbles, the total computational error accumulates. Therefore information about the error is 
crucial, since it gives us a lower bound for the total computational error. 

Two important observations illustrating the weakness of the algorithm should be noted here: 

• The density of the tessellation points near the symmetry axis {y ~ 0) increases with each time step (fig. 
[3]b). Since the time step is constant, this may lead to failure of the stability condition for the linearized 
finite difference (FD) scheme applied to the nonlinear parabolic equation ([2|) due to this algorithm. 

• The opposite effect occurs near the external boundary y = h. However, the situation here is even much 
worse. In fact, there is not enough information to reconstruct the curvature and the unit vector at a point 
Ajv lying on the boundary and the algorithm, in fact, simply eliminates it. It creates instead the point 
B'p^ along the boundary which should be the next point Bn+i (fig. Elb)). See [5] for more details. 

First we should underline here that the computations were stable (the stable steady-state regime has been 
reached) for every value of the external velocity V under consideration. Note that in our computations at high 
V , the number of tessellation points has increased from 2^ to about 220 at steady-state. As expected, the worst 
situation (in the sense of computational time) occurs for the largest value of the velocity, V = 1.560796, which 
is slightly less that the critical value Vcr = Tr/{2h) predicted by the analytical solution [l]. In numbers, it takes 
about 1 hour to complete the computation in the case V < 0.5, rising to over 48 hours for V = 1.560796 (SGI 
Altix Itanium 2, 1.5GHz). The algorithm could not reach the steady-state regime at all for V > 1.560796, and 
lost physical meaning since the interface had a branch lying outside the external boundary y = h. All this 
illustrates that the existing algorithm is well organised and works according to expectations but it is naturally 
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sensitive to the value of the boundary velocity V. Thus it makes sense to ask about algorithm accuracy for a 
specific velocity versus space and time steps. 

As the exact analytical solution to the Mullins problem is known, we can estimate errors in the computations 
for all the physical and geometrical quantities: position of the interface x{y), curvature K{y) and the velocity 
components vi{y),V2{y)- Corresponding relative errors for all solution parameters are presented in fig. |4]for 
different applied velocities. 
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Figure 4: Relative error, compared to Mullins' solution (fT7|) and (PD|) . in (a) the position of the interface x{y), 
(b) its curvature ^(y), and the velocity components (c) wi(y) and (d) V2{y), for different velocities V of the 
external boundary. 



As expected, the least accurate solutions are those for the largest velocity V ~ 1.560796. The error can be 
as large as 90% near the boundary (y = h) in the x component of velocity vi and 60% near the symmetry axis 
for the y component. The latter error is naturally related to the fact that the value of the velocity is equal to 
zero at the symmetry point. However, the drastic difference in value for this numerical noise, and its distance 
from the axis, indicates that at the larger values of V it exceeds reasonable expectations and is really related 
to the computational accuracy. 

One can also consider that the error near the external boundary is due to paramctrisation of the numerical 
solution in y rather than x, but the standard numerical algorithm tries to preserve edge lengths. Moreover, the 
algorithm introduces new points in a regular fashion. 

Thus both the errors (near the interface ends) are a consequence of the phenomena discussed in the two 
observations above. As V decreases, the accuracy increases for given bounds on the edge lengths L. 

At first glance, it would appear that the position of the interface, x{y), should be computed with better 
accuracy than all other solution parameters, which are, in fact, results of some derivative procedure. However, 
our computations show that this is not in the case and the relative error for x{y) varies from 6% to 28% for the 
velocity V = 1.560796 while the curvature error is lower. Note also that the error for smaller velocities reaches 
a few percent near the boundary or symmetry axis. 

The maximal absolute errors for all solution parameters mostly appear near the external boundary y = h{= 
1). This highlights that the implementation of the boundary condition in the existing algorithm cannot be 
considered as sufficient and should be improved. 
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Moreover, in the case of MuUins' solution an additional simple local indicator defined by identity 
(independent of the integration of the solution variables) could equally be considered. It is clear from the 
results presented in fig. Oa), that the error in this condition is not localized near the ends of the interface, as 
one might expect from the above. Moreover, this "internal error" is present for all values of V and is comparable 
with that near the interface ends. 



-0.005 



-0.01 



V=0.001 
V=0.1 
V=0.5 
V=1.0 

V=1 .560796 




0.4 0.5 0.6 
y position 



(a) 



0.006 
0.005 
0.004 


Reference case ^ 
Decrease fimestep by factor 1 x 
Haive fimestep and minimum lengfti + 

Constant lengtli 


+ 


0.003 






0.002 






0.001 


-0.001 
-0.002 






j^^»^«<»*jill»lifWMTl+.|fWlllM^^^ 

Xx 

x 




-0.003 






0.004 







0.2 0.4 0.6 O.i 

y position 



(b) 



Figure 5: Relative error for the characteristic relation (|13p for (a) varying velocity V, and (b) varying numerical 
parameters {At, L). The best accuracy is obtained by keeping the line segments of equal length. 



To investigate accurately this "internal" error we repeated the computations for a specific velocity V = 1 and 
decreased both the time increment. At, and the minimum edge length, Lmin (fig- Elb)). This has improved the 
quality of the computations, but there is still an error at some internal points of the interval that is comparable 
with the error near the ends. Unfortunately, it considerably increases the computational time by several hours. 
The obvious route of decreasing At but fixing Lmin, to ameliorate this effect, leads to greater error in the 
solution (fig. [5l[b)). The accuracy of the solution can be improved internally by making the line segments of 
equal length (fig. Eljb)) [8], although this doesn't affect the error at the boundary. 

Possible sources of the "internal" error include: a) non-optimal distribution of the tessellation points along 
the interface after some time; b) imperfections in the correction procedure (which adds and eliminates points 
from the interface at some prescribed time); c) point-to-point error variation related to the fact that the 
"diffusion" -type coefficient changes from point to point along the interface (recall that equation ([2]) is a nonlinear 
parabolic equation which is solved by a direct FD scheme with a fixed time step). 

In the last computation in figure [S]Jb) , for the line segments of equal length, we have redistributed points 
to make the segments equal at every time step, but note that this length may change in time. Apart from the 
fact that the number of tessellation points is smaller than for the standard algorithm, such a comparison is 
not absolutely fair as the redistribution in the standard algorithm was done every 20 time steps. To discover 
if there is an effect of redistribution frequency on the accuracy, we have tested these two algorithms under the 
same strategy by redistributing the points (in a different way) every time step and every 20**^ time step. The 
error in the function F defined in (|25p . which is a constant in the case of MuUins solution, are presented in 
[Gfa). It is evident that the standard algorithm is quite sensitive to the chosen strategy. For a given position 
of the points on the interface, the relative error may differ by as much as two orders of magnitude, whereas 
this is not the case for the equal segment strategy. In this case, only near the external boundary is there 
some small fluctuation in the accuracy. Comparing the two redistribution algorithms for the same frequency 
of redistribution, the largest error always appears in the case of the standard algorithm - by up to two orders 
of magnitude - despite the fact that the number of tessellation points was greater. Moreover, in the standard 
algorithm this error is irregularly distributed along the interface. Recall that the number of tessellation points 
in the standard algorithm changes during the computations from 2^ initially to around 150 (for V = 1) in the 
steady-state regime, while the number of the points in the second (equal length) algorithm remains constant. 
Therefore the computational time for the second algorithm was less by a factor of approximately two. On the 
other hand, the difference in the computational time between the different frequencies for redistribution for the 
equal length algorithm was only a few percent. This indicates the further possibility to optimise this algorithm 
by redistributing points every M time steps. It is clear that M = M(V) and this needs further investigation 
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Note that the function F from ([25|) can be used as an indicator of the accuracy of the computation only 
for the MuUins solution. However, there are three universal indicators which can be helpful to estimate the 
accuracy for any computations. Namely, the relative errors of the numerical representations of the identities 
([26)) . ((25)) and The respective results are shown in[6tb)-(d)). 
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Figure 6: Relative errors in the computations shown by the integral measures for unit velocity of the external 
boundary, V = 1.0, obtained with four different computational strategies: the standard redistribution of the 
tessellation points, tU and tl20; uniform length segment redistribution, VI and V20, at every time step and 
every 20*^ step, respectively, (a) the function F defined in ([25]) (as in Fig. [5]); (b,c,d) the three general internal 
measures; for the notation used in the graphs see Appendix IA.4I All integrals were been computed with the 
trapezium rule. 



As in the case of the specific indicator discussed above (Fig. [6]^ a)), all general indicators presented in Fig. [6] 
(b) - (d) indicate that the equidistant distribution of the tessellation points is much better than the standard 
algorithm, regardless of the chosen strategy. Moreover, even near the symmetry point, y = 0, where the value of 
the indicators all tend to zero and have a large influence on the relative errors, the accuracy of the computations 
for the second algorithm is extremely high. This is not the case for the standard algorithm. 

In Fig. [3 the relative errors of the solution variables are presented for both algorithms: the standard one 
and the equidistant distribution. The result shown in Fig. [7l[a) looks surprisingly at first glance: although the 
accuracy of the computations performed with these two algorithms is of the same order and the error related to 
the new algorithm is distributed more uniformly, it appears that the accuracy of the standard algorithm is better 
than the equal-segment-length algorithm, at least with respect to the accuracy of the position of the interface. 
However, this not in the case. In fact, as was shown above, the computational error for the standard algorithm is 
redistributed along the interface irregularly whereas that for the equal-segment algorithm is practically uniform. 
As a result, the criterion to stop the iteration process to find the steady-state solution works differently for the 
two algorithms. The prescribed maximal growth 10^^ in each time-step, measured on the axis of symmetry, is 
reached more quickly for the new algorithm. This is the second reason (together with number of tessellation 
points) why this algorithm is faster. If one were to run both algorithms for the same time, or for the same 
number of iterations, and compare the corresponding results, the discussed paradox should not appear and the 
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new algorithm always provides better accuracy by as much as two orders of magnitude. 

However, for the accuracy of other problem variables, the interface curvature, k, and the interface velocities, 
Vi and cf. in Fig. [7] (b)-(d). the new algorithm is more accurate, notwithstanding the above argument. 

Finally, let us stress again that the proposed three indicators are more versatile measures than a comparison 
of numerical steady-state solution with the analytical one, since the latter comparison includes an additional 
error related to the determination of the steady-state regime, while the indicators show us accuracy of the 
solution at any time step. 
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Figure 7: Relative error in (a) the position of the interface x{y), (b) the curvature K{y), and (c) and (d) the 
velocity components vi(y) and V2{y), for velocity of the external boundary V ~ 1 and different numerical 
algorithms. 



4 Discussion and Conclusions 

All these results clearly indicate that the existing algorithm should be used with caution, especially when 
investigating foam behaviour near the critical velocity. Moreover, when there are many bubbles in a simulation 
(fig- IHc)), the user is restricted to some critical number of tessellation points M, which gives a limitation 
on accuracy even for low velocity. In fact, the foam structure is highly non-uniform, in the sense that the 
bubbles may have different sizes. Effectively this means that every interface has its own critical velocity and 
the bigger bubbles thus introduce larger errors. This creates the following duality: to accurately describe the 
process of bubble motion it is necessary to have the computational error as small as possible, while the error 
generated near the critical velocity takes its greatest value. This is true for boundary or internal bubbles equally. 
Problems requiring high accuracy of the solution near the external boundary arc related to the investigation of 
the boundary effects describing the total phcnomcnological behaviour of the foam structures. 

Another important remark is that the choice of the initial condition for testing the numerical procedure here 
(a straight-line interface) is much more severe than any of the instantaneous solutions reported in the Appendix. 
One could even think of worse situations to test the algorithm, for example if the initial interface is not convex 
or even not smooth. This may lead to super-critical velocities and so on. 
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To revise and improve the existing numerical algorithm, we propose to use another strategy for the redistri- 
bution of the tessellation points: an equal-segment-length distribution of the tessellation points is much more 
favourable [5]. However, this strategy is not sufhcient to eliminate inaccuracy in the computation near the 
maximum velocity in the steady-state regime. The reason for this is the behaviour of the steady-state solution 
near the wall (pS)) . In fact, there exist two possible realisations of this algorithm. One, which we have used in 
our computations, is to keep the same number of tessellation points, M, then with time the length between the 
consecutive points, L, will increase significantly when V is near Vcr- This leads to an effective loss of accuracy. 
Another strategy would be to keep the same distance between points during the computations. However, the 
number of the points, M, will then increase to infinity as Vcr. Formally this should preserve computational 
accuracy but will lead to a drastic increase in computational time and memory use, which is unacceptable. Thus, 
further adaptations to the algorithm are required if high accuracy is required, for example in the steady-state 
regime with velocities near the critical value. 

Taking advantage of the auxiliary identities ([26|) . ([28|) and ([31]) . we may correct the instantaneous solution 
obtained within any algorithm at any or even every time step without time-consuming computations, as the 
identities are valid for any instantaneous solution. They also make possible further investigation of the asymp- 
totic behaviour of the bounded interface solution near the ends. For example, any possible solution behaves at 
the symmetry axis according to (|36p . which allows us to tackle the error in the solution near the symmetry axis. 
On the other hand, the results obtained in section [2731 mav allow us to construct and implement a new numerical 
procedure/elements tackling the boundary condition in a more accurate way (without losing any near-boundary 
points). 

Finally, as we have shown, the identities ((26)) . ((28)) and ([3T|) may be used to probe the accuracy of compu- 
tations. These indicators are extremely helpful as they are not based on information about the exact solution 
and can therefore illuminate inaccuracy of the numerical solution without any preliminary knowledge about the 
exact solution itself. 

Summarizing, we have shown that an improvement of the numerical algorithm is highly desirable and 
possible. Apart from the fact that some of the improvements have been indicated and proven in this paper, 
there is still an open question how to deal with accuracy of the computations near the critical velocities and 
near the external boundaries. We have also suggested possible directions for future investigation: improved 
implementation of the boundary condition, creation of additional near-boundary points. Further, to check new 
results related to the numerical algorithm we need a larger set of analytical benchmarks. 
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A Appendix: A family of symmetrical instantaneous solutions 



In this section we present analytical representations of some instantaneous symmetrical solutions for the interface 
satisfying the same boundary ([221) and symmetry (|15p conditions as MuUins' solution. 

A.l First example 

Let we consider the following simple combination of compatible velocities 

V 



v2{y) = ^W^y, v^{y) = W^\-W^y\ W ^ ^===, (39) 

which satisfy (j3ip with C = and, as a result, can be used to construct a symmetrical instantaneous solution. 
Here W is the same constant as in ((36|). Natural restrictions ([37]) for the existence of such solution give the 
same estimate: 

1 VI 
W<-, or < -, (40) 

which holds true for any values of V and h. The shape of the interface is an ellipse described by equation ([^5)1 : 

x(y, t) = x{0, t) - ^l-W^y^. (41) 

The tangential angle 9 for this solution is a monotonically decreasing function in the interval y e (0, h) and 

0{y) e {e,mn, 7r/2), e^^n = tt/2 - arcsin(iy/i) > 0. (42) 

A. 2 Second example 

Let we now consider another specific instantaneous solution assuming that vi — W < V. Then the second 
component of the velocity satisfies the equation 



To find V2{y) it is more convenient to return to the differential equation ()30p rather than working with the 
nonlinear integral equation (|43|) . After integration it takes the form 



^ (^) = -^y' (44) 



where the odd function <i> is defined as 

m = I (axctane + ^) , ^'(0 = (45) 

Note that <I> : IR+ — > [0, tt/2) is a monotonic function. Moreover, one can easily obtain the constraint Wh < 7r/4, 
which is similar to ((57)) and (|i(I)) . Then the required velocity component V2 can be found from: 

V2 = -W'i>~\Wy), (46) 

and we can finally find the complete solution using (|27p : 

Wy 

W 

Finally, note that the tangential angle 9 for this solution is a monotonically decreasing function in the 
interval y £ (0, h) 

Wh 

0{y) G (e™n,7r/2), e,mn = '^/2 - / (l + (*"'(0)') (48) 



x{y,t)^x{0,t)) + ^ I ^-\Odt (47) 
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It remains only to find possible values of the unknown constant W in order to satisfy the boundary condition 
([23|) . The relevant equation takes the form 



<^-\Wh) ^ ^[^1. (49) 

This equation has the unique solution W = W^,{V, h) < V. In fact, the left hand side is an increasing function 
from zero to infinity as Wh ^ 7r/4, whereas the right-hand side is a decreasing function taking values between 
oo when Wh and when Wh Vh. Additionally one can conclude from this that Wh < min{Fft,, 7r/4}, 
so the restriction defined after pS)) always holds. In other words, this solution, as well as that of the first 
example, is well-defined for arbitrary velocity V and position of the boundary y — h. We can also show that 
Orain IS always positive: 

7r/4 oo oo 

e™„ >^/^- J (l + = 7r/2 - y (1 + ry2) $'(^)rf^ J JT^ = ^^^'^ 



It is interesting to note that in the case « 1 both of the instantaneous solutions constructed above coincide 
with MuUins' solution to within an accuracy of 0{V^) for any fixed value of h. On the other hand, in this case 
the solution is practically (with the same accuracy) a straight line. 



A. 3 General case 

The second example above indicates how to build a wider class of symmetrical instantaneous solutions. Let us 
introduce the following set 21 C C^([— a, a]), (a > 0) of even functions, -0(0 = '0(— 0; satisfying the following 
four conditions: 

m = le + 0{e), ^^0; 0(a) <i; 0-' > 0; (^-^=^==^ > 0, ^e(0,a). (51) 

Note that a may differ from function to function, but it is necessary that for every function there exists some 
a > for which all conditions (jSip hold. For some functions it may happen that a = oo. 
For example, the following three functions belong to the set 21: 



MO = Jsin^e, MO = le, MO = I *-'(C)dC = \ii Attt^ ' 



(52) 



with ai = 7r/2, 02 = 1 and 03 = 00 respectively. These three functions have been collected from MuUins' 
solution and two previous examples. Thus, the set 2t is not empty. 

Using any function from this set we can construct a symmetrical instantaneous solution with velocity com- 
ponents in the form: 

V2[y)=-W^P'[Wy), v,{y) = Wi:'{Wy)^l^^^, (53) 

that identically satisfies equation ([?T|) with c = 0. Then the unknown constant W should be taken to be of the 
form W = W*{Vh)/h where W^,{Vh) > is a solution of the implicit equation: 

nw.) Vh ^^^^ 



v/20(VF,)(l-2V'(I^*)) 

Because of the last condition in (|5ip . there may exist only one solution of this equation. If, in addition, the 
left-hand side of ([5T|) tends to infinity as W* a, then the solution always exists and W^, < a. However, if the 
left-hand side of ((5T|) tends to a finite value i > as W» a, then the solution exists only under the additional 
condition 

Vh < La. (55) 
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One can check that for Mulhns' solution L = 1 and coincides with (fT5|) . For the other two cases previously 

discussed above, we have L = oo so the solution of the implicit equation (j54|l always exists and no solvability 

condition is needed in these cases. 

To reconstruct the complete symmetrical instantaneous solution based on (|53p it is enough to substitute it 

in dig), dm and ([Mil- 
Note that in the case V « I, the solution to gives W^, ~ V, as one can conclude from the first part of 

([5T|l . This means that any instantaneous solution differs negligibly from the Mullins' solution for small values 

of the velocity V. 

It remains to investigate two important constraints (P7|) . Taking into account that 

h h 

J V2{m = -i^iWh), J F{Od^ = i (arcsin(4V'(W^/i) - l) + f ) , 



then the two constraints ([37|) are equivalent in this case and correspond to tp{Wh) < 1/2, which coincides with 
the second part of ([5T|) . 

In fact, the third condition, -0' > 0, from (jSip was never used. We added it only to have a convex instanta- 
neous solution; without it, we can construct non-convex interfaces. 



A. 4 Integral measures 

Here we define explicitly the integral measures of problem quantities, used for assessing the accuracy of the 
computations, and their relative errors. 
After (Ull): 

intF = - J F{^)dC, relative error ^ - ij j / (^0{y) - ^tt^ - 1. (56) 

\o / 

After (EHl): 

y / y 

ii^tv= / ^r^^^' relative error = ( / ^^d^ | / {x{0) - x{y)) - I. (57) 

After (BII): 



y / y \ 

intv2 = / ..(Ode; relative error = / v.^C / " ^ (^8) 
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